Strong excitations of a Bose— Einstein condensate: Barrier resonances 



Juan J. G. Ripoll and Victor M. Perez-Garci'a 
Departamento de Matemdticas, Escuela Tecnica Superior de Ingenieros Industriales 
Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain 



oo 
on 

ON 



> 

ON 

(N 

o 
in 
o 

00 
On 

S3 
Ok 
^— > ■ 

G ■ 

s : 



We study the dynamics of a forced condensed atom cloud 
and relate the behavior to a classical Mathieu oscillator in a 
singular potential. It is found that there are wide resonances 
which can strongly affect the dynamics even when dissipation 
is present. The behavior is characteristic of condensed clouds 
of any shape and has experimental relevance. 



The recent experimental realization of Bosc-Einstein 
condensation (BEC) in ultra-cold atomic gases |l| has 
triggered the theoretical exploration of the properties of 
Bose gases. The current model used to describe a sys- 
tem with a fixed mean number N of weakly interacting 
bosons, trapped in a parabolic potential V(r) when the 
particle density and temperature of the condensate are 
small enough is the so called Gross-Pitaevskii equation 



(GPE) 



iMti) = -— V 2 V + V(r,t)iP + U a \ip\ 2 ip, 
2m 



(1) 



where Uq = 4irh 2 a/m is defined in terms of the ground 
state scattering length a. The normalization for ip 
is N = J \rb\ 2 d 3 r, and the trapping potential is 
given by V(r,t) = \mv 2 (\ 2 x (t)x 2 + \ 2 {t)y 2 + \ 2 z (t)z 2 ). 
X v , (j) — x, y, z) are, as usual, functions that describe the 
anisotropics of the trap Q . In real experiments with sta- 
tionary systems they are constants and the geometry of 
the trap imposes the condition \ x = X y = 1. X z = v : z / 'v 
is the quotient between the trap frequency along the z- 
direction v z and the radial one v r = v. 

The problem of the excitations of the condensate under 
a periodic driving was studied experimentally in Ref . || . 
On the theoretical side an analysis in the Thomas-Fermi 
limit which providing an approximation to the low-energy 
excitation spectra were done in Ref. Q and numerical 
simulations done in Ref. [^). More accurate theoretical 
predictions were found using time dependent variational 
methods Q. Other results have been found using dif- 
ferent approaches to the problem for the single and 
double condensate Q cases. Recent research on this area 
has focused also in the effects of damping on the spec- 
trum of low energy excitations PJT(| . 

Most of the theoretical work done up to now is devoted 
to the analysis of the condensate properties for weak per- 
turbations. It is our intention here to extend the analysis 
to the case when the perturbation is stronger or is main- 
tained for a longer time than in the first experiments || . 



An unexpected result is the existence of a parametric res- 
onance which has important experimental implications. 

Let us start by deriving our model equations. It can 
be proved that every solution of Eq. (Q) is a stationary 
point of an action corresponding, up to a divergence, to 
the Lagrangian density C: 

£ = y - rdti>) 

+^-\Vi/j\ 2 + V(rM 2 + U^\\ (2) 

where the asterisk denotes complex conjugation. That 
is, instead of working with the GPE we can treat the 
action, S = J CcPrdt = J t f L(t)dt, and study its in- 
variance properties and extrema. For instance, from the 
invariance of (|2|) under global phase transformations, one 
can assure the conservation the number of particles in the 
Bose condensed state N = J \ip\ 2 d 3 r. 

To further simplify the problem, we restrict the shape 
of the function ip to a convenient family of trial functions 
and study the time evolution of the parameters that de- 
fine that family. A natural choice, as discussed in is 
a three dimensional Gaussian-like function 



$(T,t)=A n 



-lv- 



lo] 




+ir;a, +ir) /3, ( 



(3) 



We use this trial function to obtain an averaged La- 
grangian and then the Lagrange equations to find the 
evolution equations for all the parameters. Of special rel- 
evance are the equations for the widths, which after defin- 
ing the constants P — y // 2/irNa/ao and ag — ^Jfij ' [mv), 
as well as a set of rescaled variables for time, r = z/t, and 
the widths, = aov Vl {rj = x, y, z), are found to be 



1 



v x + \ x {t)v x = — -f- 



P 



y 



P 



,, M 1 P 
v z + \ 2 z (t)v z = — + 2. 



(4a) 
(4b) 
(4c) 



As can be seen, the variational equations remain the same 
as those of || with the only change of time-dependent A 
values. Although a lot of different dynamics are possible 
in Eqs. ( pajjic] ), we will restrict ourselves to the param- 
eter ranges relevant for current Bose-Einstein condensa- 
tion experiments. 
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To gain insight on the problem let us first analyze the 
radially symmetric version of Eq. (0), which is 



1 



+ \\t)v = — + 



P 

7a' 



(5) 



This equation is a Hill's equation with a singular poten- 
tial given by the right hand side of (||). Approaching 
the experimental setups, we will treat the case where the 
potential strength is harmonic, i.e. X 2 (t) = 1 + ecos(uit). 

First, as it can be seen in Fig. ^, if one reproduces 
the evolution of the Gaussian condensate with a far from 
equilibrium width in a stationary trap (e = 0), what 
comes out is a periodic orbit with a frequency that is 
slightly different from the natural one, and with fast 
bounces near the origin. In other words, for large ampli- 
tude motion, the origin acts as an elastical barrier and the 
harmonic frequencies of the trap gain importance over 
the details of the potential well and its linearization. 

Let us now concentrate on the time dependent prob- 
lem (e ^ 0). We have solved it numerically, as stated in 
Eq. ([5]). Scanning the parameter space (u, e) one finds at 
least two wide resonance regions where the radial width 
v grows exponentially as far as the perturbation is main- 
tained (Fig. ||(a)). Both regions have the form of wedges, 
with a starting point, (u! m in, £min), and a growing width 
as e is increased (Fig. ||(b)). 

Exhaustive simulations have also been made for sev- 
eral values of the parameter P going from P = 9.2 -the 
JILA H| experiment- to 20 times this value, and replac- 
ing the singular potential with others of a similar shape 
(l/i> 4 , l/i> 3 , etc). The description is the same as the one 
stated in the preceding paragraphs, the base frequencies 
-u m in = 2 and Lu rn in = 1- remaining the same up to a 
0.5% precision. 

A great care is needed when treating Eq. (||) numeri- 
cally because of the singularity at v = 0. We have used 
a Dormand-Prince pair |TJ], the ODE Suite of MAT- 
LAB, and Vazquez's conservative scheme H|. When one 
approaches a resonance condition all those methods fail 
after a sufficiently long run and we have employed sev- 
eral stiff methods: the BDF formulas @], the LSODE 
Fortran library and the stiff integrators included in the 
MATLAB ODE suite. 

Since the origin acts as an elastic wall it is intuitively 
appealing to replace the singular potential in (|J) with a 
bounce condition on the origin, i.e. an impact oscillator 
p| . We replace Eq. (||) with the following one 



lim (u, v) 

t—>ta 



+ \ 2 {t)v 
(0 + ,K) 



0, 

lim (v, v) 
t->tt 



(6) 



- (c\+ 



(o + ,-v c ), 



where t c denotes any isolated instance when the system 
bounces against the v = singularity. Let us show that 
this equation is in turn equivalent to an elastic oscillator 
without barrier conditions. We introduce the change of 



variables v — \u\, where u is an unrestricted real function 
which satisfies the following one-dimensional harmonic 
oscillator equation 



+ X 2 (t)u = 0, 



(7) 



This equation is thus a Mathieu equation which arises 
in the analysis of parametrically forced oscillators. It is 
a well know problem where one can analytically obtain 
jfjj a lot of information. It is now easy to prove that 
every solution of Eq. (0) provides a solution of Eq. (|^). 
And vice versa, from every solution of Eq. (|^) it is pos- 
sible to construct a solution of Eq. (^), unique up to a 
sign. Among the properties of the Mathieu equation the 
one that we are concerned most about is the existence 
of instability regions in the parameter space. The lim- 
its of these zones can be found by means of Floquet's 
theory, and have the shape of wedges that start on the 
points (u m i n , e min ) = (2, 0), (1, 0), (2/3, 0), and widen 
as e is increased up from zero, in close similarity to our 
numerical results from Eq. (g). 

Either with an asymptotic method, or by making use 
of the singular perturbation theory, it is possible to study 
the evolution of the condensate around the resonances. 
For a perturbation frequency close enough to the first 
resionance, that is for|w — 2 1 = 2<5 = o(l), an asymptotic 
method fl2|| yields, up to first order, 



± 



4^2 



<5 2 , 



u(t) ~ ce qt cos(ut + e ). 



(8a) 
(8b) 



Here we see that for some values of S and e the expo- 
nent q is a positive real number and the amplitude of the 
oscillations grows unlimitedly. Also, the strength of the 
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resonance is maximum for a value of uj max ~ 2 — ^- and a 
second order Taylor expansion in Eq. ( |Sa|) allows to find 
the approximate resonance conditions \u> — 2| < | + 
which is in agreement with the numerical simulations of 
(||) . The treatment of other resonances is more difficult 
as they are caused by higher order terms; this means 
that they have a smaller region of influence and are not 
so strong. We have observed at least two subharmonic 
resonances in our numerical simulations. 

The resonances here found resist the presence of dis- 
sipation. Adding a viscous damping term to Eq. (^) it 
becomes 



u + (1 + e cos(ojt))u + 7« = 0. 



(9) 



As before, it is also possible to find the approximate form 
of the dominant contribution on resonance 



u(t) ~ ce (9 - 7)t cos(^ + 6» ), 



(10) 



where q is given by an equation similar to (|8a|). Due 
to 7, the resonance regions in the parameter space are 
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constrained to values of (u>, e) for which the strength of 
the resonance, q, is greater that of the dissipative term. 

For instance, taking the data from the JILA experi- 
ment H , we can estimate a condensate lifetime of about 
110 ms and a value for 7 of 0.15 in natural units of the 
condensate. Such a damping makes the e TO j„ value raise 
from 0.09 to 0.18 for the P = 9.2 case. Thus, the insta- 
bility should not be appreciated unless the perturbation 
amplitude exceeds the 20%. 

An interesting effect of damping is that the evolution 
of a continuously perturbed condensate outside the in- 
stability regions develops a simple limit cycle perfectly 
synchronized to the frequency of the parametric pertur- 
bation, and with a size that depends only on the pertur- 
bation parameters, (w, e). The optimal frequency for the 
creation of the limit cycle depends on the amplitude of 
the perturbation. For very small perturbations the fre- 
quency is close to that of the small amplitude oscillations. 
However, as the perturbation is increased, the frequency 
rapidly approaches the Mathieu resonance. The appear- 
ance of a limit cycle opens the door to a wide family 
of phenomena, from chaotic motion to bifurcation the- 
ory. This limit cycle would exist under a great variety of 
dissipative terms, and is not exclusive of linear damping. 

Apart from studying the variational reduction of the 
radially symmetric problem, we have also studied numer- 
ically the exact radially symmetric version ofEq. § In 
a rather complete inspection of the parameter space, we 
have found that, up to the point in which the numerical 
simulations could be continuated, the resonances exist 
and behave like the variational model predicted. This 
is a somewhat surprising result since the partial differ- 
ential equation has a lot of degree of freedom available 
and one could think that the resonant gaussian-like mode 
would decay into a combination of nonresonant modes. 
Although some contribution of higher order modes is gen- 
erated by the perturbation the energy pumping mech- 
anism provided by the resonance is quite efficient and 
seemingly active as far as the perturbation is maintained. 
More details of the analysis will be provided elsewhere. 

Finally let us comment on the nonsymmetric case ruled 
by the full set of Eqs. ( ^aj^c| ). Following the experimen- 
tal setups P], we can once more take a sinusoidal time 
dependence for every A^(t) coefficient: 

^(*) = + e„ cos(urf)), r) = x,y, z. (11) 

This choice accounts both for the m = (e x = e y , e z = 0) 
and the m = 2 (e x = —e y , e z = 0) perturbations from the 
JILA experiment Q . In the latter case the potential is a 
parabolic one, with fixed frequencies on a rotating frame. 

Substituting our effective perturbation frequencies into 
Eq. (^) we get a set of three coupled Mathieu equations 
with a potential that is singular on the v x = 0, v y = and 
v z — planes. The singularities are at least as strong as 
1/v 3 , and the numerical simulations again confirm that 



they act as elastical walls, so we now proceed with the 
change of variables = to find 

u v +\2(t)u ri =Q (12) 

for 77 = x, y, z. Now the situation is a bit more complex. 
The first new feature we find is the existence of several 
sets of instability regions. Due to having three a priori 
different constants Ao^, the three oscillators in Eq. jl^ ) 
are not equivalent and we may get three sets of reso- 
nances in the (e^ , uj) space. Further details will be given 
in a future publication jlj] . 

A very interesting result can be found by analyzing 
exactly the center of mass movement. Defining < 77 >^= 
J rj\ifj\ 2 dr, and computing its time derivatives using Eq. 
(0) we arrive again to a Mathieu equation for the center 
of mass Q 

d 2 1 

^2 < V >i> +2"^^) < V >i>= (13) 

with rj = x,y, z. This has a serious consequence which is 
that feeding the condensate will result in the exponential 
amplification of any initial displacement of the center of 
mass. This is an exact prediction based only on the GPE. 
On the other hand, we can imagine two ways in which this 
effect could have experimental relevance. First, the study 
of the strength of the resonance should account for possi- 
ble dissipation effects due to collisions. And secondly, it 
would be interesting to study how the condensed and the 
noncondensed centers of mass response to the perturba- 
tion, because their different dissipation regimes may lead 
to an effective separation of both clouds. It also imposes 
restrictions on the time the perturbation can be applied 
if the condensate is to remain stable. 

Summarizing, we have found that medium to large 
amplitude oscillations of a condensate approach the har- 
monic trap frequencies, not the ones resulting from the 
linearization of the variational equations. Even when 
damping is added only frequencies close to the Mathieu 
resonance regions do excite the condensate as a whole in 
an efficient way, causing the appearance of a stable limit 
cycle. 

We have also found that for this kind of parametrical 
drive the resonances are naturally wide in all perturba- 
tion regimes for this kind of parametrical drive. This 
width grows with the strength of the interaction, a fact 
that can be checked in the experiments by forcing the 
system for a longer time than what it is currently done. 

The variational method showed that the kinetic terms 
in the evolution equations guarantee a 1/v 3 singularity 
as far as we impose a repulsive interaction between the 
atoms in the cloud. This singularity is enough to cause 
the appearance of the instability regions. Thus, we have 
a wide family of systems that behave much the same. 
Also the response of the noncondensed atoms under the 
parametrical perturbation will be qualitatively similar to 
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that of the condensed ones, with the only difference that 
the former are subject to a more intense dissipation. This 
dissipation can be enough to distinguish both kind of flu- 
ids: while the condensed part might suffer an exponential 
growth, the uncondensed part might develop bounded os- 
cillations. 

This work has been supported in part by the Spanish 
Ministry of Education and Culture under grants PB95- 
0389, PB96-0534 and AP97-08930807. 



FIG. 2. Example of resonance (a) Time dependent behav- 
ior of v for P = 9.2, lo = 2.04 and e = 0.15 and initial 
conditions (y,v) = (1.6,0) close to the equilibrium point (b) 
One of the regions of instability (shaded) in the (lo, e) plane 
corresponding to the resonance lo = 2.04 

FIG. 3. Plot of the evolution of a cylindrically symmet- 
ric condensate for P = 9.2 under a sinusoidal perturbation 
(lo, e) = (2.04,0.1) of the radial strength of the trap. Both 
(a) the axial v z and (b) the radial v r widths are plotted. 
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FIG. 1. Orbits in the phase space evolution for the vari- 
ational approximation of a spherically symmetric condensate 
with P = 9.2, e = 0, and several different initial conditions. 
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